{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# point_defect_diffusion calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "Description updated: 2019-07-26\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The point_defect_diffusion calculation style estimates the diffusion rate of a point defect at a specified temperature.  A system is created with a single point defect, and subjected to a long time molecular dynamics simulation.  The mean square displacement for the defect is computed, and used to estimate the diffusion constant.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- The calculation estimates the defect's diffusion by computing the mean square displacement of all atoms in the system.  This is useful for estimating rates associated with vacancies and self-interstitials as the process is not confined to a single atom's motion.  However, this makes the calculation ill-suited to measuring diffusion of substitutional impurities as it does not individually track each atom's position throughout.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "First, a defect system is constructed by adding a single point defect (or defect cluster) to an initially bulk system using the atomman.defect.point() function.\n",
    "\n",
    "A LAMMPS simulation is then performed on the defect system.  The simulation consists of two separate runs\n",
    "\n",
    "1. NVT equilibrium run: The system is allowed to equilibrate at the given temperature using nvt integration.\n",
    "\n",
    "2. NVE measurement run: The system is then evolved using nve integration, and the total mean square displacement of all atoms is measured as a function of time.\n",
    "\n",
    "Between the two runs, the atomic velocities are scaled such that the average temperature of the nve run is closer to the target temperature.\n",
    "\n",
    "The mean square displacement of the defect, $\\left< \\Delta r_{ptd}^2 \\right>$ is then estimated using the mean square displacement of the atoms $\\left< \\Delta r_{i}^2 \\right>$.  Under the assumption that all diffusion is associated with the single point defect, the defect's mean square displacement can be taken as the summed square displacement of the atoms\n",
    "\n",
    "$$ \\left< \\Delta r_{ptd}^2 \\right> \\approx \\sum_i^N \\Delta r_{i}^2 = N \\left< \\Delta r_{i}^2 \\right>, $$\n",
    "\n",
    "where $N$ is the number of atoms in the system.  The diffusion constant is then estimated by linearly fitting the change in mean square displacement with time\n",
    "\n",
    "$$ \\left< \\Delta r_{ptd}^2 \\right> = 2 d D_{ptd} \\Delta t, $$\n",
    "\n",
    "where d is the number of dimensions included.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [DataModelDict](https://github.com/usnistgov/DataModelDict)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2019-07-29 using iprPy version 0.9.0\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import sys\n",
    "import uuid\n",
    "import shutil\n",
    "import random\n",
    "import datetime\n",
    "from copy import deepcopy\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np \n",
    "\n",
    "# https://github.com/usnistgov/DataModelDict \n",
    "from DataModelDict import DataModelDict as DM\n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'point_defect_diffusion'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ (required) is the LAMMPS command to use.\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential_file__ gives the path to the potential_LAMMPS reference record to use.  Here, this parameter is automatically generated using potential_name and librarydir.\n",
    "\n",
    "- __potential_dir__ gives the path for the folder containing the artifacts associated with the potential (i.e. eam.alloy file).  Here, this parameter is automatically generated using potential_name and librarydir.\n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required).  Here, this parameter is automatically generated from potential_file and potential_dir."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully loaded potential 1999--Mishin-Y--Ni--LAMMPS--ipr1\n"
     ]
    }
   ],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Define potential_file and potential_dir using librarydir and potential_name\n",
    "potential_file = Path(iprPy.libdir, 'potential_LAMMPS', f'{potential_name}.json')\n",
    "potential_dir = Path(iprPy.libdir, 'potential_LAMMPS', potential_name)\n",
    "\n",
    "# Initialize Potential object using potential_file and potential_dir.\n",
    "potential = lmp.Potential(potential_file, potential_dir)\n",
    "print('Successfully loaded potential', potential)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __prototype_name__ gives the name of the crystal_prototype reference record in the iprPy library to load. \n",
    "\n",
    "- __symbols__ is a list of the potential's elemental model symbols to associate with the unique atom types of the loaded system. \n",
    "\n",
    "- __box_parameters__ is a list of the a, b, c lattice constants to assign to the loaded file.\n",
    "\n",
    "- __load_file__ gives the path to the atomic configuration file to load for the ucell system.  Here, this is generated automatically using prototype_name and librarydir.\n",
    "\n",
    "- __load_style__ specifies the format of load_file.  Here, this is automatically set for crystal_prototype records.\n",
    "\n",
    "- __load_options__ specifies any other keyword options for properly loading the load_file.  Here, this is automatically set for crystal_prototype records.\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is generated using the load parameters and symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.520,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.520,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.520]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.760 |   1.760\n",
      "      2 |       1 |   1.760 |   0.000 |   1.760\n",
      "      3 |       1 |   1.760 |   1.760 |   0.000\n"
     ]
    }
   ],
   "source": [
    "prototype_name = 'A1--Cu--fcc'\n",
    "symbols = ['Ni']\n",
    "box_parameters = uc.set_in_units([3.52, 3.52, 3.52], 'angstrom')\n",
    "\n",
    "# Define load_file using librarydir and prototype_name\n",
    "load_file = Path(iprPy.libdir, 'crystal_prototype', f'{prototype_name}.json')\n",
    "\n",
    "# Define load_style and load_options for crystal_prototype records\n",
    "load_style = 'system_model'\n",
    "load_options = {}\n",
    "\n",
    "# Create ucell by loading prototype record\n",
    "ucell = am.load(load_style, load_file, symbols=symbols, **load_options)\n",
    "\n",
    "# Rescale ucell using box_parameters\n",
    "ucell.box_set(a=box_parameters[0], b=box_parameters[1], c=box_parameters[2], scale=True)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Specify the defect parameters\n",
    "\n",
    "- __pointdefect_name__ gives the name of a point-defect reference record in the iprPy library containing point defect input parameters.\n",
    "\n",
    "- __pointdefect_file__ gives the path to a point_defect reference containing point defect input parameters.  Here, this is built automatically using pointdefect_name and librarydir.\n",
    "\n",
    "- __point_kwargs__ (required) is a dictionary or list of dictonaries containing parameters for generating the defect. Here, values are extracted from pointdefect_file. Allowed keywords are:\n",
    "\n",
    "    - __ptd_type__ indicates which defect type to generate: 'v' for vacancy, 'i' for interstitial, 's' for substitutional, or 'db' for dumbbell.\n",
    "    \n",
    "    - __atype__ is the atom type to assign to the defect atom ('i', 's', 'db' ptd_types).\n",
    "    \n",
    "    - __pos__ specifies the position for adding the defect atom (all ptd_types).\n",
    "    \n",
    "    - __ptd_id__ specifies the id of an atom in the initial system where the defect is to be added. Alternative to using pos ('v', 's', 'db' ptd_types).\n",
    "    \n",
    "    - __db_vect__ gives the vector associated with the dumbbell interstitial to generate ('db' ptd_type).\n",
    "    \n",
    "    - __scale__ indicates if pos and db_vect are in absolute (False) or box-relative (True) coordinates. Default is False.\n",
    "    \n",
    "    - __atol__ is the absolute tolerance for position-based searching. Default is 1e-3 angstroms.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "point_kwargs =\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[{'ptd_type': 'v', 'pos': array([0., 0., 0.]), 'scale': False}]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pointdefect_name = 'A1--Cu--fcc--vacancy'\n",
    "#pointdefect_name = 'A1--Cu--fcc--1nn-divacancy'\n",
    "#pointdefect_name = 'A1--Cu--fcc--2nn-divacancy'\n",
    "#pointdefect_name = 'A1--Cu--fcc--100-dumbbell'\n",
    "#pointdefect_name = 'A1--Cu--fcc--110-dumbbell'\n",
    "#pointdefect_name = 'A1--Cu--fcc--111-dumbbell'\n",
    "#pointdefect_name = 'A1--Cu--fcc--octahedral-interstitial'\n",
    "#pointdefect_name = 'A1--Cu--fcc--tetrahedral-interstitial'\n",
    "#pointdefect_name = 'A1--Cu--fcc--crowdion-interstitial'\n",
    "\n",
    "# Define pointdefect_file using librarydir and pointdefect_name\n",
    "pointdefect_file = Path(iprPy.libdir, 'point_defect', f'{pointdefect_name}.json')\n",
    "\n",
    "# Parse pointdefect_file using iprPy.input.interpret()\n",
    "defectinputs = {'ucell':ucell, 'pointdefect_file':pointdefect_file}\n",
    "iprPy.input.subset('pointdefect').interpret(defectinputs)\n",
    "\n",
    "# Extract point_kwargs\n",
    "point_kwargs = defectinputs['point_kwargs']\n",
    "print('point_kwargs =')\n",
    "point_kwargs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Modify system\n",
    "\n",
    "- __sizemults__ list of three integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in creating system.\n",
    "\n",
    "- __system__ is an atomman.System to perform the scan on (required). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# of atoms in system = 4000\n"
     ]
    }
   ],
   "source": [
    "sizemults = [10, 10, 10]\n",
    "\n",
    "# Generate system by supersizing ucell\n",
    "system = ucell.supersize(*sizemults)\n",
    "print('# of atoms in system =', system.natoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.6. Specify calculation-specific run parameters\n",
    "\n",
    "- __temperature__ temperature in Kelvin at which to run the MD integration scheme at.  Default value is '0'.\n",
    "\n",
    "- __runsteps__ specifies how many timesteps to integrate the system.  Default value is 200000.\n",
    "\n",
    "- __thermosteps__ specifies how often LAMMPS prints the system-wide thermo data.  Default value is runsteps/1000, or 1 if runsteps is less than 1000.\n",
    "    \n",
    "- __dumpsteps__ specifies how often LAMMPS saves the atomic configuration to a LAMMPS dump file.  Default value is runsteps, meaning only the first and last states are saved.\n",
    "    \n",
    "- __equilsteps__ specifies how many timesteps are ignored as equilibration time when computing the mean box parameters.  Default value is 10000.\n",
    "\n",
    "- __randomseed__ provides a random number seed to generating the initial atomic velocities.  Default value gives a random number as the seed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = 1900\n",
    "runsteps = 200000\n",
    "thermosteps = 100\n",
    "dumpsteps = runsteps\n",
    "equilsteps = 20000\n",
    "randomseed = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. diffusion.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('diffusion.template', 'w') as f:\n",
    "    f.write(\"\"\"# LAMMPS input script for dynamic msd computation\n",
    "\n",
    "box tilt large\n",
    "\n",
    "<atomman_system_info>\n",
    "\n",
    "<atomman_pair_info>\n",
    "\n",
    "# Assign simulation parameter values\n",
    "variable temperature equal <temperature>\n",
    "variable randomseed equal <randomseed>\n",
    "variable thermosteps equal <thermosteps>\n",
    "variable timestep equal <timestep>\n",
    "variable equilsteps equal <equilsteps>\n",
    "variable dumpsteps equal <dumpsteps>\n",
    "variable runsteps equal <runsteps>\n",
    "variable twotemp equal 2*${temperature}\n",
    "variable damptemp equal 100*${timestep}\n",
    "\n",
    "# Specify property computes\n",
    "compute peatom all pe/atom\n",
    "compute msd all msd com yes\n",
    "\n",
    "# Define thermo data\n",
    "thermo ${thermosteps}\n",
    "thermo_style custom step temp pe pxx pyy pzz c_msd[1] c_msd[2] c_msd[3] c_msd[4]\n",
    "thermo_modify format float %.13e\n",
    "\n",
    "# Specify timestep\n",
    "timestep ${timestep}\n",
    "\n",
    "# Create velocities and equilibrate system using nvt\n",
    "velocity all create ${twotemp} ${randomseed}\n",
    "fix 1 all nvt temp ${temperature} ${temperature} ${damptemp}\n",
    "run ${equilsteps}\n",
    "unfix 1\n",
    "<dump_info>\n",
    "\n",
    "# Scale velocities to wanted temperature and run nve\n",
    "velocity all scale ${temperature}\n",
    "reset_timestep 0\n",
    "fix 2 all nve\n",
    "run ${runsteps}\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2. pointdiffusion()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pointdiffusion(lammps_command, system, potential, point_kwargs,\n",
    "                   mpi_command=None, temperature=300,\n",
    "                   runsteps=200000, thermosteps=None, dumpsteps=0,\n",
    "                   equilsteps=20000, randomseed=None):\n",
    "                   \n",
    "    \"\"\"\n",
    "    Evaluates the diffusion rate of a point defect at a given temperature. This\n",
    "    method will run two simulations: an NVT run at the specified temperature to \n",
    "    equilibrate the system, then an NVE run to measure the defect's diffusion \n",
    "    rate. The diffusion rate is evaluated using the mean squared displacement of\n",
    "    all atoms in the system, and using the assumption that diffusion is only due\n",
    "    to the added defect(s).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    point_kwargs : dict or list of dict\n",
    "        One or more dictionaries containing the keyword arguments for\n",
    "        the atomman.defect.point() function to generate specific point\n",
    "        defect configuration(s).\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    temperature : float, optional\n",
    "        The temperature to run at (default is 300.0).\n",
    "    runsteps : int, optional\n",
    "        The number of integration steps to perform (default is 200000).\n",
    "    thermosteps : int, optional\n",
    "        Thermo values will be reported every this many steps (default is\n",
    "        100).\n",
    "    dumpsteps : int or None, optional\n",
    "        Dump files will be saved every this many steps (default is 0,\n",
    "        which does not output dump files).\n",
    "    equilsteps : int, optional\n",
    "        The number of timesteps at the beginning of the simulation to\n",
    "        exclude when computing average values (default is 20000).\n",
    "    randomseed : int or None, optional\n",
    "        Random number seed used by LAMMPS in creating velocities and with\n",
    "        the Langevin thermostat.  (Default is None which will select a\n",
    "        random int between 1 and 900000000.)\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'natoms'** (*int*) - The number of atoms in the system.\n",
    "        - **'temp'** (*float*) - The mean measured temperature.\n",
    "        - **'pxx'** (*float*) - The mean measured normal xx pressure.\n",
    "        - **'pyy'** (*float*) - The mean measured normal yy pressure.\n",
    "        - **'pzz'** (*float*) - The mean measured normal zz pressure.\n",
    "        - **'Epot'** (*numpy.array*) - The mean measured total potential \n",
    "          energy.\n",
    "        - **'temp_std'** (*float*) - The standard deviation in the measured\n",
    "          temperature values.\n",
    "        - **'pxx_std'** (*float*) - The standard deviation in the measured\n",
    "          normal xx pressure values.\n",
    "        - **'pyy_std'** (*float*) - The standard deviation in the measured\n",
    "          normal yy pressure values.\n",
    "        - **'pzz_std'** (*float*) - The standard deviation in the measured\n",
    "          normal zz pressure values.\n",
    "        - **'Epot_std'** (*float*) - The standard deviation in the measured\n",
    "          total potential energy values.\n",
    "        - **'dx'** (*float*) - The computed diffusion constant along the \n",
    "          x-direction.\n",
    "        - **'dy'** (*float*) - The computed diffusion constant along the \n",
    "          y-direction.\n",
    "        - **'dz'** (*float*) - The computed diffusion constant along the \n",
    "          y-direction.\n",
    "        - **'d'** (*float*) - The total computed diffusion constant.\n",
    "    \"\"\"\n",
    "    try:\n",
    "        # Get script's location if __file__ exists\n",
    "        script_dir = Path(__file__).parent\n",
    "    except:\n",
    "        # Use cwd otherwise\n",
    "        script_dir = Path.cwd()\n",
    "\n",
    "    # Add defect(s) to the initially perfect system\n",
    "    if not isinstance(point_kwargs, (list, tuple)):\n",
    "        point_kwargs = [point_kwargs]\n",
    "    for pkwargs in point_kwargs:\n",
    "        system = am.defect.point(system, **pkwargs)\n",
    "    \n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "    \n",
    "    #Get lammps version date\n",
    "    lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "    \n",
    "    # Check that temperature is greater than zero\n",
    "    if temperature <= 0.0:\n",
    "        raise ValueError('Temperature must be greater than zero')\n",
    "    \n",
    "    # Handle default values\n",
    "    if thermosteps is None: \n",
    "        thermosteps = runsteps // 1000\n",
    "        if thermosteps == 0:\n",
    "            thermosteps = 1\n",
    "    if dumpsteps is None:\n",
    "        dumpsteps = runsteps\n",
    "    if randomseed is None:\n",
    "        randomseed = random.randint(1, 900000000)\n",
    "    \n",
    "    # Define lammps variables\n",
    "    lammps_variables = {}\n",
    "    system_info = system.dump('atom_data', f='initial.dat',\n",
    "                              units=potential.units,\n",
    "                              atom_style=potential.atom_style)\n",
    "    lammps_variables['atomman_system_info'] = system_info\n",
    "    lammps_variables['atomman_pair_info'] = potential.pair_info(system.symbols)\n",
    "    lammps_variables['temperature'] = temperature\n",
    "    lammps_variables['runsteps'] = runsteps\n",
    "    lammps_variables['equilsteps'] = equilsteps\n",
    "    lammps_variables['thermosteps'] = thermosteps\n",
    "    lammps_variables['dumpsteps'] = dumpsteps\n",
    "    lammps_variables['randomseed'] = randomseed\n",
    "    lammps_variables['timestep'] = lmp.style.timestep(potential.units)\n",
    "    \n",
    "    # Set dump_info\n",
    "    if dumpsteps == 0:\n",
    "        lammps_variables['dump_info'] = ''\n",
    "    else:\n",
    "        lammps_variables['dump_info'] = '\\n'.join([\n",
    "            '',\n",
    "            '# Define dump files',\n",
    "            'dump dumpit all custom ${dumpsteps} *.dump id type x y z c_peatom',\n",
    "            'dump_modify dumpit format <dump_modify_format>',\n",
    "            '',\n",
    "        ])\n",
    "        \n",
    "        # Set dump_modify_format based on lammps_date\n",
    "        if lammps_date < datetime.date(2016, 8, 3):\n",
    "            lammps_variables['dump_modify_format'] = '\"%d %d %.13e %.13e %.13e %.13e\"'\n",
    "        else:\n",
    "            lammps_variables['dump_modify_format'] = 'float %.13e'\n",
    "    \n",
    "    # Write lammps input script\n",
    "    template_file = Path(script_dir, 'diffusion.template')\n",
    "    lammps_script = 'diffusion.in'\n",
    "    with open(template_file) as f:\n",
    "        template = f.read()\n",
    "    with open(lammps_script, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables, '<', '>'))\n",
    "    \n",
    "    # Run lammps\n",
    "    output = lmp.run(lammps_command, 'diffusion.in', mpi_command)\n",
    "    \n",
    "    # Extract LAMMPS thermo data.\n",
    "    thermo = output.simulations[1]['thermo']\n",
    "    temps = thermo.Temp.values\n",
    "    pxxs = uc.set_in_units(thermo.Pxx.values, lammps_units['pressure'])\n",
    "    pyys = uc.set_in_units(thermo.Pyy.values, lammps_units['pressure'])\n",
    "    pzzs = uc.set_in_units(thermo.Pzz.values, lammps_units['pressure'])\n",
    "    potengs = uc.set_in_units(thermo.PotEng.values, lammps_units['energy'])\n",
    "    steps = thermo.Step.values\n",
    "    \n",
    "    # Read user-defined thermo data\n",
    "    if output.lammps_date < datetime.date(2016, 8, 1):\n",
    "        msd_x = uc.set_in_units(thermo['msd[1]'].values,\n",
    "                                lammps_units['length']+'^2')\n",
    "        msd_y = uc.set_in_units(thermo['msd[2]'].values,\n",
    "                                lammps_units['length']+'^2')\n",
    "        msd_z = uc.set_in_units(thermo['msd[3]'].values,\n",
    "                                lammps_units['length']+'^2')\n",
    "        msd = uc.set_in_units(thermo['msd[4]'].values,\n",
    "                              lammps_units['length']+'^2')\n",
    "    else:\n",
    "        msd_x = uc.set_in_units(thermo['c_msd[1]'].values,\n",
    "                                lammps_units['length']+'^2')\n",
    "        msd_y = uc.set_in_units(thermo['c_msd[2]'].values,\n",
    "                                lammps_units['length']+'^2')\n",
    "        msd_z = uc.set_in_units(thermo['c_msd[3]'].values,\n",
    "                                lammps_units['length']+'^2')\n",
    "        msd = uc.set_in_units(thermo['c_msd[4]'].values,\n",
    "                              lammps_units['length']+'^2')\n",
    "        \n",
    "    # Initialize results dict\n",
    "    results = {}\n",
    "    results['natoms'] = system.natoms\n",
    "    \n",
    "    # Get mean and std for temperature, pressure, and potential energy\n",
    "    results['temp'] = np.mean(temps)\n",
    "    results['temp_std'] = np.std(temps)\n",
    "    results['pxx'] = np.mean(pxxs)\n",
    "    results['pxx_std'] = np.std(pxxs)\n",
    "    results['pyy'] = np.mean(pyys)\n",
    "    results['pyy_std'] = np.std(pyys)\n",
    "    results['pzz'] = np.mean(pzzs)\n",
    "    results['pzz_std'] = np.std(pzzs)\n",
    "    results['Epot'] = np.mean(potengs)\n",
    "    results['Epot_std'] = np.std(potengs)\n",
    "    \n",
    "    # Convert steps to times\n",
    "    times = steps * uc.set_in_units(lammps_variables['timestep'], lammps_units['time'])\n",
    "    \n",
    "    # Estimate diffusion rates\n",
    "    # MSD_ptd = natoms * MSD_atoms (if one defect in system)\n",
    "    # MSD = 2 * ndim * D * t  -->  D = MSD/t / (2 * ndim)\n",
    "    mx = np.polyfit(times, system.natoms * msd_x, 1)[0]\n",
    "    my = np.polyfit(times, system.natoms * msd_y, 1)[0]\n",
    "    mz = np.polyfit(times, system.natoms * msd_z, 1)[0]\n",
    "    m = np.polyfit(times, system.natoms * msd, 1)[0]\n",
    "    \n",
    "    results['dx'] = mx / 2\n",
    "    results['dy'] = my / 2\n",
    "    results['dz'] = mz / 2\n",
    "    results['d'] = m / 6\n",
    "    \n",
    "    return results\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1. Generate point defect system and evaluate the energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = pointdiffusion(lammps_command, system, potential, point_kwargs,\n",
    "                              mpi_command = mpi_command,\n",
    "                              temperature = temperature,\n",
    "                              runsteps = runsteps,\n",
    "                              thermosteps = thermosteps,\n",
    "                              dumpsteps = dumpsteps,\n",
    "                              equilsteps = equilsteps,\n",
    "                              randomseed = randomseed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['natoms', 'temp', 'temp_std', 'pxx', 'pxx_std', 'pyy', 'pyy_std', 'pzz', 'pzz_std', 'Epot', 'Epot_std', 'dx', 'dy', 'dz', 'd'])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1 Define units for outputting values\n",
    "\n",
    "- __length2_pertime_unit__ is the units to display the computed diffusion constants in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "length2_pertime_unit = 'm^2/s'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2 Print directional diffusion rates, $D_x$, $D_y$, and $D_z$ and total diffusion rate, $D$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dx = 1.289386392173063e-10 m^2/s\n",
      "Dy = 1.8253069489194243e-10 m^2/s\n",
      "Dz = 1.2757396287269496e-10 m^2/s\n",
      "D =  1.4634776566064903e-10 m^2/s\n"
     ]
    }
   ],
   "source": [
    "print('Dx =', uc.get_in_units(results_dict['dx'], length2_pertime_unit), length2_pertime_unit)\n",
    "print('Dy =', uc.get_in_units(results_dict['dy'], length2_pertime_unit), length2_pertime_unit)\n",
    "print('Dz =', uc.get_in_units(results_dict['dz'], length2_pertime_unit), length2_pertime_unit)\n",
    "print('D = ', uc.get_in_units(results_dict['d'], length2_pertime_unit), length2_pertime_unit)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
